library(knitr)
# purl("08_fec_multivariate.Rmd", output = "08_fec_multivariate.R", documentation = 2)
rerun_MCMC <- FALSE
set.seed(224819)
library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0 ✔ purrr 0.2.5
## ✔ tibble 1.4.2 ✔ dplyr 0.7.6
## ✔ tidyr 0.8.1 ✔ stringr 1.3.1
## ✔ readr 1.1.1 ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(parallel)
source("color_map.R")
source("ggplot_theme.R")
h2fec <- read.table("../../Data/Processed/eggs_per_vial.txt",
sep = '\t',
dec = ".", header = TRUE,
stringsAsFactors = FALSE)
# Standardize egg_total
h2fec$egg_total <- as.numeric(scale(h2fec$egg_total))
h2fec$animal <- seq(1, nrow(h2fec))
h2fec$treat <- as.factor(h2fec$treat)
pedigree <- h2fec[, c("animal", "sireid", "damid")]
names(pedigree) <- c("animal", "sire", "dam")
pedigree$animal <- as.character(pedigree$animal)
pedigree$sire <- as.character(pedigree$sire)
pedigree$dam <- as.character(pedigree$dam)
sires <- data.frame(animal = unique(pedigree$sire),
sire = NA, dam = NA, stringsAsFactors = FALSE)
dams <- data.frame(animal = unique(pedigree$dam),
sire = NA, dam = NA, stringsAsFactors = FALSE)
pedigree <- bind_rows(sires, dams, pedigree) %>%
as.data.frame()
genet_corr <- tibble(model = character(),
HS_STD = numeric(),
LY_STD = numeric(),
HS_LY = numeric(),
n_eff = numeric())
iter <- 6.5e6
burnin <- 5e4
thinning <- 500
chains <- 12
if (rerun_MCMC) {
HS <- h2fec %>%
filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>%
as.data.frame()
DR <- h2fec %>%
filter(treat == "LY") %>% rename(Eggs_DR = egg_total) %>%
as.data.frame()
STD <- h2fec %>%
filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>%
as.data.frame()
h2fec_mrg <- full_join(full_join(HS, DR), STD)
prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
priors <- list(prior1)
prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_DR) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2fec_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("../../Data/Processed/fec_total_multivariate_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_DR:traitEggs_DR.animal"])
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
DR_STD = median(DR_STD),
HS_DR = median(HS_DR),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "../../Data/Processed/Genetic_Correlations_Fecundity.csv")
}
load("../../Data/Processed/fec_total_multivariate_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)
plot(fe[, 2, drop = FALSE], ask = FALSE)
plot(fe[, 3, drop = FALSE], ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)
plot(re[, 2, drop = FALSE], ask = FALSE)
plot(re[, 3, drop = FALSE], ask = FALSE)
autocorr.diag(re)
## traitEggs_STD:traitEggs_STD.animal
## Lag 0 1.0000000000
## Lag 500 0.1874099388
## Lag 2500 0.0368637247
## Lag 5000 0.0068407312
## Lag 25000 -0.0009086811
## traitEggs_HS:traitEggs_STD.animal
## Lag 0 1.0000000000
## Lag 500 0.3135652226
## Lag 2500 0.0488681915
## Lag 5000 0.0004422247
## Lag 25000 0.0066357983
## traitEggs_DR:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 0.245970917
## Lag 2500 0.044196514
## Lag 5000 0.011865920
## Lag 25000 0.002057849
## traitEggs_STD:traitEggs_HS.animal
## Lag 0 1.0000000000
## Lag 500 0.3135652226
## Lag 2500 0.0488681915
## Lag 5000 0.0004422247
## Lag 25000 0.0066357983
## traitEggs_HS:traitEggs_HS.animal
## Lag 0 1.0000000000
## Lag 500 0.2983162660
## Lag 2500 0.0387993500
## Lag 5000 0.0007746715
## Lag 25000 -0.0006467404
## traitEggs_DR:traitEggs_HS.animal
## Lag 0 1.000000000
## Lag 500 0.329916759
## Lag 2500 0.051896081
## Lag 5000 0.002264740
## Lag 25000 0.004759115
## traitEggs_STD:traitEggs_DR.animal
## Lag 0 1.000000000
## Lag 500 0.245970917
## Lag 2500 0.044196514
## Lag 5000 0.011865920
## Lag 25000 0.002057849
## traitEggs_HS:traitEggs_DR.animal
## Lag 0 1.000000000
## Lag 500 0.329916759
## Lag 2500 0.051896081
## Lag 5000 0.002264740
## Lag 25000 0.004759115
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.343334829 0.1694615436
## Lag 2500 0.071567781 0.0343858227
## Lag 5000 0.012262294 0.0066624941
## Lag 25000 -0.003595793 -0.0006694499
## traitEggs_HS.units traitEggs_DR.units
## Lag 0 1.0000000000 1.0000000000
## Lag 500 0.1354468964 0.1653264427
## Lag 2500 0.0181847424 0.0320639066
## Lag 5000 -0.0004421368 0.0057609407
## Lag 25000 -0.0045816742 0.0007743094
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## 74878.16 59041.49
## traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## 64199.26 59041.49
## traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## 64269.67 57406.62
## traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## 64199.26 57406.62
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## 50166.62 78571.11
## traitEggs_HS.units traitEggs_DR.units
## 96791.83 78128.78
# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,] 0.2730985 0.12469556
## [2,] 0.2995462 0.07412571
## [3,] 0.4452922 0.18096616
## [4,] 0.1252125 0.07514900
## [5,] 0.1048659 0.04575423
## [6,] 0.1438339 0.05751000
## [7,] 0.2007193 0.03386499
## traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,] 0.24849777 0.12469556
## [2,] 0.20893779 0.07412571
## [3,] 0.20795884 0.18096616
## [4,] 0.11505847 0.07514900
## [5,] 0.08302472 0.04575423
## [6,] 0.08359697 0.05751000
## [7,] 0.13166049 0.03386499
## traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## [1,] 0.110130511 0.11774301
## [2,] 0.034842647 0.04079747
## [3,] 0.118196674 0.07136862
## [4,] 0.070105046 0.08114866
## [5,] 0.022681420 0.03957909
## [6,] 0.035602734 0.03912530
## [7,] 0.008878417 0.02152333
## traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## [1,] 0.24849777 0.11774301
## [2,] 0.20893779 0.04079747
## [3,] 0.20795884 0.07136862
## [4,] 0.11505847 0.08114866
## [5,] 0.08302472 0.03957909
## [6,] 0.08359697 0.03912530
## [7,] 0.13166049 0.02152333
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## [1,] 0.23588453 0.4891190
## [2,] 0.16338372 0.4610870
## [3,] 0.10745186 0.4486223
## [4,] 0.11301484 0.6035127
## [5,] 0.07074388 0.6815414
## [6,] 0.05895181 0.4694439
## [7,] 0.09380421 0.5364391
## traitEggs_HS.units traitEggs_DR.units
## [1,] 0.1884776 0.2815059
## [2,] 0.1863606 0.2473274
## [3,] 0.1944373 0.3105254
## [4,] 0.1768785 0.3693298
## [5,] 0.2189036 0.3670518
## [6,] 0.2347433 0.4200357
## [7,] 0.2539427 0.3827891
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])
plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)
median(HS_STD)
## [1] 0.7124564
HPDinterval(HS_STD)
## lower upper
## var1 -0.5112175 0.9987634
## attr(,"Probability")
## [1] 0.95
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_DR:traitEggs_DR.animal"])
plot(DR_STD)
median(DR_STD)
## [1] 0.9312949
HPDinterval(DR_STD)
## lower upper
## var1 0.3639826 0.99934
## attr(,"Probability")
## [1] 0.95
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_DR)
median(HS_DR)
## [1] 0.7238309
HPDinterval(HS_DR)
## lower upper
## var1 -0.5517535 0.9979309
## attr(,"Probability")
## [1] 0.95
save(re, file = "../../Data/Processed/fec_total_multivariate_model_output.Rda")
library(tidyverse)
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
##
## ggsave
B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`DR vs. STD` = as.numeric(DR_STD),
`HS vs. DR` = as.numeric(HS_DR))
B <- as_tibble( B ) %>%
select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>%
rename(`HS vs. C` = `HS vs. STD`) %>%
rename(`DR vs. C` = `DR vs. STD`) %>%
rename(`HS vs. DR` = `HS vs. DR`)
colMeans(B)
## HS vs. C DR vs. C HS vs. DR
## 0.5413373 0.8384703 0.5400363
B %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 1) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.10, 0.85),
text = element_text(size = 10),
legend.text = element_text(size = 10))+
scale_x_continuous(expand = c(0, 0), limits = c(-1.2, 1.2)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 12)) +
my_theme
ggsave(last_plot(), file = "../../Figures/Fecundity_Genetic_Correlation_Plot.pdf",
width = 4, height = 4)
Fecundity_correlation <- last_plot()
save(Fecundity_correlation, file = "../../Figures/Fecundity_correlation.Rda")
if (rerun_MCMC) {
iter <- 1e5
burnin <- 2e3
thinning <- 500
chains <- 12
HS <- h2fec %>%
filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>%
as.data.frame()
DR <- h2fec %>%
filter(treat == "LY") %>% rename(Eggs_DR = egg_total) %>%
as.data.frame()
STD <- h2fec %>%
filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>%
as.data.frame()
h2fec_mrg <- full_join(full_join(HS, DR), STD)
prior1 <- list(R = list(V = diag(3), nu = 1.002, fix = 3),
G = list(G1 = list(V = diag(3),
nu = 1,
alpha.mu = rep(0, 3),
alpha.V = diag(3) * 25)))
priors <- list(prior1)
prior_names <- c("Tri: Parameter expanded")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_DR) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2fec_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("../../Data/Processed/fec_total_multivariate_model_prior_expanded", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_DR:traitEggs_DR.animal"])
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
DR_STD = median(DR_STD),
HS_DR = median(HS_DR),
n_eff = mean(n_eff)))
}
write_csv(
genet_corr,
path = "../../Data/Processed/Genetic_Correlations_Fecundity_Parameter_Expanded.csv")
}
load("../../Data/Processed/fec_total_multivariate_model_prior_expanded1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)
plot(fe[, 2, drop = FALSE], ask = FALSE)
plot(fe[, 3, drop = FALSE], ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)
plot(re[, 2, drop = FALSE], ask = FALSE)
plot(re[, 3, drop = FALSE], ask = FALSE)
autocorr.diag(re)
## traitEggs_STD:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 -0.003725192
## Lag 2500 0.019490364
## Lag 5000 0.030377401
## Lag 25000 -0.001895760
## traitEggs_HS:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 -0.016525153
## Lag 2500 -0.013641051
## Lag 5000 0.002471646
## Lag 25000 0.010064237
## traitEggs_DR:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 0.028389065
## Lag 2500 0.006133787
## Lag 5000 -0.040036034
## Lag 25000 -0.001075109
## traitEggs_STD:traitEggs_HS.animal
## Lag 0 1.000000000
## Lag 500 -0.016525153
## Lag 2500 -0.013641051
## Lag 5000 0.002471646
## Lag 25000 0.010064237
## traitEggs_HS:traitEggs_HS.animal
## Lag 0 1.0000000000
## Lag 500 0.0403265957
## Lag 2500 0.0024874149
## Lag 5000 -0.0045582535
## Lag 25000 0.0003642594
## traitEggs_DR:traitEggs_HS.animal
## Lag 0 1.0000000000
## Lag 500 -0.0108965197
## Lag 2500 0.0198383220
## Lag 5000 -0.0001302109
## Lag 25000 0.0080958169
## traitEggs_STD:traitEggs_DR.animal
## Lag 0 1.000000000
## Lag 500 0.028389065
## Lag 2500 0.006133787
## Lag 5000 -0.040036034
## Lag 25000 -0.001075109
## traitEggs_HS:traitEggs_DR.animal
## Lag 0 1.0000000000
## Lag 500 -0.0108965197
## Lag 2500 0.0198383220
## Lag 5000 -0.0001302109
## Lag 25000 0.0080958169
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## Lag 0 1.000000000 1.000000000
## Lag 500 0.023406187 -0.022469930
## Lag 2500 0.010900801 0.013481221
## Lag 5000 -0.004557295 0.020502363
## Lag 25000 0.011769764 -0.006948712
## traitEggs_HS.units traitEggs_DR.units
## Lag 0 1.000000000 NaN
## Lag 500 -0.005884207 NaN
## Lag 2500 0.008485809 NaN
## Lag 5000 0.014211774 NaN
## Lag 25000 -0.010895854 NaN
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## 2437.161 2464.272
## traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## 2381.271 2464.272
## traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## 2444.515 2431.644
## traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## 2381.271 2431.644
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## 2303.084 2460.700
## traitEggs_HS.units traitEggs_DR.units
## 2372.235 0.000
# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,] 0.3590307 0.1282773450
## [2,] 0.2442366 0.0006156665
## [3,] 0.1861954 -0.0100121837
## [4,] 0.4102305 0.1024781219
## [5,] 0.1929804 0.0546066988
## [6,] 0.2738059 0.0193737220
## [7,] 0.2613179 0.0017872135
## traitEggs_DR:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,] -0.0083137086 0.1282773450
## [2,] 0.1011918879 0.0006156665
## [3,] 0.0018253463 -0.0100121837
## [4,] 0.0001729217 0.1024781219
## [5,] 0.0564094448 0.0546066988
## [6,] -0.0117165610 0.0193737220
## [7,] 0.0178162082 0.0017872135
## traitEggs_HS:traitEggs_HS.animal traitEggs_DR:traitEggs_HS.animal
## [1,] 7.439954e-02 -0.0051252069
## [2,] 1.955109e-05 0.0005158809
## [3,] 6.917970e-02 -0.0003176388
## [4,] 5.379830e-02 0.0036724221
## [5,] 8.522121e-02 0.0503100316
## [6,] 4.261395e-02 0.0069020223
## [7,] 6.829332e-02 0.0013074032
## traitEggs_STD:traitEggs_DR.animal traitEggs_HS:traitEggs_DR.animal
## [1,] -0.0083137086 -0.0051252069
## [2,] 0.1011918879 0.0005158809
## [3,] 0.0018253463 -0.0003176388
## [4,] 0.0001729217 0.0036724221
## [5,] 0.0564094448 0.0503100316
## [6,] -0.0117165610 0.0069020223
## [7,] 0.0178162082 0.0013074032
## traitEggs_DR:traitEggs_DR.animal traitEggs_STD.units
## [1,] 0.004483110 0.4588685
## [2,] 0.058540298 0.3533715
## [3,] 0.000445515 0.5218645
## [4,] 0.001503260 0.1974659
## [5,] 0.046957532 0.6654494
## [6,] 0.003940158 0.4780136
## [7,] 0.001747412 0.4619696
## traitEggs_HS.units traitEggs_DR.units
## [1,] 0.2339418 1
## [2,] 0.3095291 1
## [3,] 0.1936399 1
## [4,] 0.2018433 1
## [5,] 0.1806931 1
## [6,] 0.1818133 1
## [7,] 0.1943753 1
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])
plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)
median(HS_STD)
## [1] 0.3946399
HPDinterval(HS_STD)
## lower upper
## var1 -0.7472588 0.9999947
## attr(,"Probability")
## [1] 0.9498299
# STD vs. DR
DR_STD <- re[ , "traitEggs_DR:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_DR:traitEggs_DR.animal"])
plot(DR_STD)
median(DR_STD)
## [1] 0.405907
HPDinterval(DR_STD)
## lower upper
## var1 -0.8619617 0.9999912
## attr(,"Probability")
## [1] 0.9498299
# DR vs. HS
HS_DR <- re[ , "traitEggs_HS:traitEggs_DR.animal"] /
sqrt(re[ , "traitEggs_DR:traitEggs_DR.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_DR)
median(HS_DR)
## [1] 0.2382352
HPDinterval(HS_DR)
## lower upper
## var1 -0.9068489 0.9999722
## attr(,"Probability")
## [1] 0.9498299
library(tidyverse)
library(cowplot)
B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`DR vs. STD` = as.numeric(DR_STD),
`HS vs. DR` = as.numeric(HS_DR))
B <- as_tibble( B ) %>%
select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>%
rename(`HS vs. C` = `HS vs. STD`) %>%
rename(`DR vs. C` = `DR vs. STD`) %>%
rename(`HS vs. DR` = `HS vs. DR`)
colMeans(B)
## HS vs. C DR vs. C HS vs. DR
## 0.2945234 0.2643091 0.1457471
B %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 1) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.10, 0.85),
text = element_text(size = 10),
legend.text = element_text(size = 10))+
scale_x_continuous(expand = c(0, 0), limits = c(-1, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 12)) +
my_theme
ggsave(last_plot(), file = "../../Figures/Fecundity_Genetic_Correlation_Plot.pdf",
width = 4, height = 4)
Fecundity_correlation <- last_plot()
save(Fecundity_correlation, file = "../../Figures/Fecundity_correlation.Rda")